* QCEW_models_xsec.do
* 2016.06.23
* Last update 2017.06.27: update path locals
* Estimates cross-sectional reduced-form models using QCEW wage bills

capture log close
set more off
timer clear 1
timer on 1
clear
set matsize 10000
set scheme s1mono

local work "/DIRECTORY"

log using "`work'/logs/QCEW_models_xsec.log", replace
eststo drop *


*****************************
* Flow control
*****************************
local descriptive = 1 /* descriptive tables */
local models = 1 /* linear models */
local robustness = 1
local graphs = 1 /* reduced-form type graphs */
local output = 1


*****************************
* Startup
*****************************
* Reading data
use "`work'/data/qcew/QCEW_xsec.dta", clear

* Subsetting to QCEW, continental US
keep if inlist(tzcode, 2, 3, 5, 6)

* Variable handling
mkspline lat_ 10 = latitude
mkspline pd_ 5 = pop_density
mkspline cd_ 3 = coast_dist
gen medage2 = median_age^2
label variable sunset_time_avg "Avg. sunset time"

* Locals for control variables
* Note only have 6 occupation shares in these data
local geographic = "coast_dist lat_1-lat_10"
local demographic = "pct_female pct_black pct_white median_age pd_*" /* educ_* pct_asian occ_* */
gen tz6Xcd = 6.tzcode#c.coast_dist
gen tz3Xcd = 3.tzcode#c.coast_dist
local controls "2.tzcode 3.tzcode 5.tzcode 6.tzcode tz3Xcd tz6Xcd" /* i.tzcode  */
local controls2 "2.tzcode 3.tzcode 5.tzcode 6.tzcode tz3Xcd tz6Xcd lat_* median_age pct_female pct_black pct_white" /* i.tzcode  */
local controls3 "2.tzcode 3.tzcode 5.tzcode 6.tzcode tz3Xcd tz6Xcd lat_* coast_dist coast cc pd_* median_age medage2 pct_female pct_black pct_white"
local educ "pct_lessthan9thgr pct_gr9to12 pct_HSgrad pct_some_college pct_BA_BS"

* Defining "donuts" around TZ boundaries
local donutPM "-117.38, -114.74"
local donutMC "-104.58, -100.58"
local donutCE "-89.28, -84.93"
local donutPM2 "-117.5, -115"
local donutMC2 "-104.5, -101"
local donutCE2 "-89.5, -85"


*****************************
* QCEW models
*****************************
if `descriptive'==1 {

* Descriptive stats
estpost summarize AverageWeeklyWage10 AverageWeeklyWage101 AverageWeeklyWage102 sunset_time, listwise
esttab using "`work'/tables/QCEW_descriptive.tex", label cells("mean sd min max") ///
	substitute(mean Mean sd Stdev min Min max Max) ///
	title("Descriptive Statistics: QCEW") ///
	nomtitle nonumber tex fragment replace
	
* Average wage by quarter - for appendix table
eststo drop *
bysort quarter: eststo: estpost summarize ln_avwage10, listwise
esttab using "`work'/tables/QCEW_wage_byqtr.tex", label cells("mean sd") ///
	substitute(mean Mean sd Stdev ln_avwage10 "Log avg. wage" " 1&" Q1& " 2&" Q2& " 3&" Q3& " 4&" Q4&) ///
	title("Mean Log Wage by Quarter: QCEW") ///
	/*nomtitle*/ nonumber nodepvar tex replace

} /* end flow control for descriptive */


if `models'==1 {
	eststo clear
	eststo: regress ln_avwage10 `controls' sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo: regress ln_avwage10 `controls2' sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo: regress ln_avwage10 `controls' 2.tzcode#c.sunset_time_avg 3.tzcode#c.sunset_time_avg 5.tzcode#c.sunset_time_avg 6.tzcode#c.sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo: regress ln_avwage10 `controls2' 2.tzcode#c.sunset_time_avg 3.tzcode#c.sunset_time_avg 5.tzcode#c.sunset_time_avg 6.tzcode#c.sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	* Table output
	esttab using "`work'/tables/QCEW_RF.tex", keep(*sunset_time_avg) order(sunset_time_avg 6.tzcode#c.sunset_time_avg 5.tzcode#c.sunset_time_avg 2.tzcode#c.sunset_time_avg 3.tzcode#c.sunset_time_avg) ///
		indicate("TZ dummies=*.tzcode" "Coastal distance=*cd" "Demographic controls=median_age") ///
		star(* 0.10 ** 0.05 *** 0.01) ///
		coeflabels(sunset_time_avg "Sunset time" 2.tzcode#c.sunset_time_avg "Central*sunset time" 3.tzcode#c.sunset_time_avg "Eastern*sunset time" 5.tzcode#c.sunset_time_avg "Mountain*sunset time" 6.tzcode#c.sunset_time_avg "Pacific*sunset time") ///
		mlabels("Log avg. wage" "Log avg. wage" "Log avg. wage" "Log avg. wage") ///
		title("Reduced-form effects on log wages") se nonumbers label tex fragment nogaps replace /*r2*/ 
} /* End flow control for `models' */


if `robustness'==1 {
	* Base spec
	eststo clear
	regress ln_avwage10 `controls3' sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_base.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Geo controls only
	eststo clear
	regress ln_avwage10 2.tzcode 3.tzcode 5.tzcode 6.tzcode tz3Xcd tz6Xcd lat_* coast_dist coast cc sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_geoonly.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		keep(sunset_time_avg) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Longitude control
	eststo clear
	regress ln_avwage10 /*2.tzcode 3.tzcode 5.tzcode 6.tzcode*/ tz3Xcd tz6Xcd lat_* coast_dist coast cc pd_* median_age pct_female pct_black pct_white longitude sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_longitude.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		keep(sunset_time_avg) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* State clustering
	eststo clear
	regress ln_avwage10 `controls3' sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(cluster statecode)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_statecluster.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Top 100 counties by population
	preserve
	eststo clear
	keep if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
	gsort -AvgEmployment10
	gen top100=(_n<=100)
	regress ln_avwage10 `controls3' sunset_time_avg if top100==1, vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_top100.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	restore
	* No eastern TZ
	eststo clear
	regress ln_avwage10 `controls3' sunset_time_avg if tzcode!=3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_noE.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Exclude San Francisco, Los Angeles, Chicago, Boston, and New York MSAs
	eststo clear
	gen SF = (fips=="06001" | fips=="06013" | fips=="06041" | fips=="06055" | fips=="06075" | fips=="06081" | fips=="06085" | fips=="06095" | fips=="06097" | fips=="06077" | fips=="06087" | fips=="06069")
	gen LA = (fips=="06037" | fips=="06059" | fips=="06111" | fips=="06071" | fips=="06065")
	gen CHI = (fips=="55059" | fips=="17097" | fips=="17111" | fips=="17037" | fips=="17089" | fips=="17043" | fips=="17031" | fips=="17197" | fips=="17093" | fips=="17063" | fips=="17091" | fips=="18089" | fips=="18111" | fips=="18073" | fips=="18127" | fips=="18091")
	gen BOS = (fips=="25025" | fips=="25021" | fips=="25023" | fips=="25017" | fips=="25009" | fips=="25027" | fips=="25005" | fips=="25001" | fips=="09015" | fips=="33015" | fips=="33017" | fips=="33011" | fips=="33013" | fips=="33001" | substr(fips, 1, 2)=="44")
	gen NYC = (fips=="36061" | fips=="36005" | fips=="36081" | fips=="36047" | fips=="36085" | fips=="36119" | fips=="36087" | fips=="36071" | fips=="36103" | fips=="36059" | fips=="36079" | fips=="36027" | fips=="36111" | fips=="09001" | fips=="09009" | fips=="09005" | fips=="42103" | fips=="42025" | fips=="42077" | fips=="42095" | fips=="42089" | fips=="34003" | fips=="34017" | fips=="34023" | fips=="34025" | fips=="34029" | fips=="34031" | fips=="34013" | fips=="34039" | fips=="34027" | fips=="34035" | fips=="34037" | fips=="34019" | fips=="34041" | fips=="34021")
	generate excludeflag=(SF==1 | LA==1 | CHI==1 | BOS==1 | NYC==1)
	regress ln_avwage10 `controls3' sunset_time_avg if excludeflag!=1 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_noHWcities.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Albouy QOL control
	replace fips2 = substr(fips, 1, 2) + "000"
	merge m:1 fips2 using "`work'/data/qol/QOL_index.dta", keep(1 3)
	eststo clear
	regress ln_avwage10 `controls3' qol sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_QOL.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' _cons qol) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
	* Educational attainment shares
	eststo clear
	regress ln_avwage10 `controls3' `educ' sunset_time_avg if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), vce(robust)
	eststo rf
	if `output' == 1 { 
		esttab rf using "`work'/tables/QCEW_linear_robust_edushares.tex", ///
		star(* 0.10 ** 0.05 *** 0.01) ///  
		drop(`controls3' `educ' _cons) ///
		se(a2) b(a2) noconstant nonumbers noobs nogaps wide label tex begin(" & ") fragment replace nomtitles nolines
	}
} /* End flow control for `robustness' */

if `graphs'==1 {
* Small changes in donuts lead to changes in default smoothing parameter, so need to hard-code smoothing parameter
local bw "1.3"
local bw2 ".16"

* Residualizing variables for plots
regress ln_avwage10 `controls'
predict rwage if e(sample), residuals
* NB below model yield statistically significant relationship (prob from non-square US)
regress sunset_time_avg `controls'
predict rSST if e(sample), residuals
regress longitude `controls'
predict rlong if e(sample), residuals

* Residualizing without overlapping longitudes
regress ln_avwage10 `controls' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rwageNO if e(sample), residuals
regress sunset_time_avg `controls' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rSSTNO if e(sample), residuals
regress longitude `controls' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rlongNO if e(sample), residuals

* Residualizing, no overlap, richer controls
regress ln_avwage10 `controls2' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rwageNO2 if e(sample), residuals
regress sunset_time_avg `controls2' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rSSTNO2 if e(sample), residuals
regress longitude `controls2' if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')
predict rlongNO2 if e(sample), residuals

* Additional variables
gen rounded_long = (-1)*round(longitude)

* Labeling variables to get better axis titles
label variable ln_avwage10 "ln(Avg. wkly wage)"
label variable rwage "ln(Avg. wkly wage), residualized"
label variable rwageNO "ln(Avg. wkly wage), residualized"
label variable rwageNO2 "ln(Avg. wkly wage), residualized"
label variable sunset_time_avg "Avg. sunset time"
label variable rSST "Avg. sunset time, residualized"
label variable rSSTNO "Avg. sunset time, residualized"
label variable rSSTNO2 "Avg. sunset time, residualized"
label variable longitude "Longitude"
label variable rlong "Longitude, residualized"
label variable rlongNO "Longitude, residualized"
label variable rlongNO2 "Longitude, residualized"

* Lpoly: log wage - longitude, separate lpoly for each TZ
twoway (lpolyci ln_avwage10 longitude if tzcode==6, degree(1) bwidth(1)) || ///
	(lpolyci ln_avwage10 longitude if tzcode==5, degree(1) bwidth(1)) || ///
	(lpolyci ln_avwage10 longitude if tzcode==2, degree(1) bwidth(1)) || ///
	(lpolyci ln_avwage10 longitude if tzcode==3, degree(1) bwidth(1)), ///
	title("") xtitle("Longitude") ytitle("ln(Avg. wkly wage)") text(5.9 -121 "Pacific" 5.9 -110 "Mountain" 5.9 -95 "Central" 5.9 -77.5 "Eastern") legend(off)
graph export "`work'/graphs/lpoly_wage_long.eps", as(eps) replace
* Binscatter: log wage - longitude
binscatter ln_avwage10 longitude if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), by(tzcode) xq(rounded_long) colors(black) ///
	legend(off) text(6.35 -119 "Pacific" 6.35 -110 "Mountain" 6.35 -95 "Central" 6.35 -77.5 "Eastern") title("") xtitle("Longitude") ytitle("ln(Avg. wkly wage)")
graph export "`work'/graphs/binscatter_wage_long.eps", as(eps) replace
* Lpoly: resid log wage - longitude
twoway (lpolyci rwage longitude if tzcode==6, degree(1)) ||  ///
	(lpolyci rwage longitude if tzcode==5, degree(1)) || ///
	(lpolyci rwage longitude if tzcode==2, degree(1)) || ///
	(lpolyci rwage longitude if tzcode==3, degree(1)), ///
	title("") xtitle("Longitude") ytitle("ln(Avg. wkly wage), residualized") text(-.4 -121 "Pacific" -.4 -110 "Mountain" -.4 -95 "Central" -.4 -77.5 "Eastern") legend(off)
graph export "`work'/graphs/lpoly_rwage_long.eps", as(eps) replace
* Lpoly: resid log wage - longitude, excluding border TZs
twoway (lpolyci rwageNO longitude if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO longitude if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO longitude if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO longitude if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')), ///
	title("") xtitle("Longitude") ytitle("ln(Avg. wkly wage), residualized") text(-.4 -121 "Pacific" -.4 -110 "Mountain" -.4 -95 "Central" -.4 -77.5 "Eastern") legend(off)
graph export "`work'/graphs/lpoly_rwage_long_nooverlap.eps", as(eps) replace
* Lpoly: resid log wage - longitude, excluding border TZs
twoway (lpolyci rwageNO2 longitude if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO2 longitude if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO2 longitude if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')) || ///
	(lpolyci rwageNO2 longitude if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw')), ///
	title("") xtitle("Longitude") ytitle("ln(Avg. wkly wage), residualized") text(-.4 -121 "Pacific" -.4 -110 "Mountain" -.4 -95 "Central" -.4 -77.5 "Eastern") legend(off)
graph export "`work'/graphs/lpoly_rwage_long_nooverlap_mc.eps", as(eps) replace
bysort tzcode: summ longitude if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), detail

* Binscatter: resid log wage - longitude, excluding border TZs
binscatter rwage longitude if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), by(tzcode) nquantiles(100) ///
	legend(off) text(.3 -120 "Pacific" .3 -110 "Mountain" .3 -95 "Central" .3 -77.5 "Eastern") xtitle("Longitude") ytitle("ln(Avg. wkly wage), residualized") ///
	msymbols(P P P P) mcolors(black black black black) lcolors(black black black black)
graph export "`work'/graphs/binscatter_rwage_long_nooverlap.eps", replace
* xq() option allows me to balance dots more evenly over TZs
binscatter rwage longitude if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), by(tzcode) xq(rounded_long) ///
	legend(off) text(.25 -120 "Pacific" .25 -110 "Mountain" .25 -95 "Central" .25 -77.5 "Eastern") xtitle("Longitude") ytitle("ln(Avg. wkly wage), residualized") ///
	msymbols(P P P P) mcolors(black black black black) lcolors(black black black black)
graph export "`work'/graphs/binscatter_rwage_long_nooverlap2.eps", replace
* Lpoly: resid log wage - resid long, no overlap
* This graph is qualitatively similar, but residualized longitude is way too hard to interpret
lpoly rwageNO rlongNO if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw') noscatter ///
	legend(off) title("") ///
	addplot(lpoly rwageNO longitude if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw') || ///
	lpoly rwageNO longitude if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) bwidth(`bw') || ///
	lpoly rwageNO longitude if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1)) bwidth(`bw')
graph export "`work'/graphs/lpoly_rwage_rlong_nooverlap.eps", as(eps) replace


* Lpoly: log wage - SST
lpoly ln_avwage10 sunset_time if tzcode==6, degree(1) noscatter legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly ln_avwage10 sunset_time if tzcode==5, degree(1) lpattern(dash) lcolor(black) || ///
	lpoly ln_avwage10 sunset_time if tzcode==2, degree(1) lpattern(dash_dot) lcolor(black) || ///
	lpoly ln_avwage10 sunset_time if tzcode==3, degree(1) lpattern(shortdash) lcolor(black))
graph export "`work'/graphs/lpoly_wage_SST.eps", as(eps) replace
* Lpoly: resid log wage - SST
lpoly rwage sunset_time if tzcode==6, degree(1) noscatter legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwage sunset_time if tzcode==5, degree(1) lpattern(dash) lcolor(black) || ///
	lpoly rwage sunset_time if tzcode==2, degree(1) lpattern(dash_dot) lcolor(black) || ///
	lpoly rwage sunset_time if tzcode==3, degree(1) lpattern(shortdash) lcolor(black))
graph export "`work'/graphs/lpoly_rwage_SST.eps", as(eps) replace
* Lpoly: resid log wage - SST, no overlap
lpoly rwageNO sunset_time if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) noscatter ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwageNO sunset_time if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(dash) lcolor(black) || ///
	lpoly rwageNO sunset_time if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(dash_dot) lcolor(black) || ///
	lpoly rwageNO sunset_time if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(shortdash) lcolor(black))
graph export "`work'/graphs/lpoly_rwage_SST_nooverlap.eps", as(eps) replace
* Lpoly: resid log wage - SST, no overlap, more controls
lpoly rwageNO2 sunset_time if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) noscatter ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwageNO2 sunset_time if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(dash) lcolor(black) || ///
	lpoly rwageNO2 sunset_time if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(dash_dot) lcolor(black) || ///
	lpoly rwageNO2 sunset_time if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), degree(1) lpattern(shortdash) lcolor(black))
graph export "`work'/graphs/lpoly_rwage_SST_nooverlap_mc.eps", as(eps) replace
* Binscatter: resid log wage - SST, no overlap, more controls
binscatter rwageNO2 sunset_time if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), by(tzcode) colors(black) ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") xtitle("Avg. sunset time") ytitle("ln(Avg. wkly wage), residualized")
graph export "`work'/graphs/binscatter_rwage_SST_nooverlap_mc.eps", as(eps) replace
* Lpoly: resid log wage - resid SST
lpoly rwage rSST if tzcode==6, degree(1) noscatter legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwage rSST if tzcode==5, degree(1) lpattern(dash) lcolor(black) || ///
	lpoly rwage rSST if tzcode==2, degree(1) lpattern(dash_dot) lcolor(black) || ///
	lpoly rwage rSST if tzcode==3, degree(1) lpattern(shortdash) lcolor(black))
graph export "`work'/graphs/lpoly_rwage_rSST.eps", as(eps) replace
* Lpoly: resid log wage - resid SST, no overlap
lpoly rwageNO rSSTNO if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') noscatter ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwageNO rSSTNO if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(dash) lcolor(black) || ///
	lpoly rwageNO rSSTNO if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') &  rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(dash_dot) lcolor(black) || ///
	lpoly rwage rSSTNO if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') &  rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(shortdash) lcolor(black)) 
graph export "`work'/graphs/lpoly_rwage_rSST_nooverlap.pdf", as(pdf) replace
* Binscatter: resid log wage - resid SST, no overlap
binscatter rwageNO rSSTNO if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE'), by(tzcode) colors(black) ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") xtitle("Avg. sunset time, residualized") ytitle("ln(Avg. wkly wage), residualized")
graph export "`work'/graphs/binscatter_rwage_rSST_nooverlap.eps", as(eps) replace
* Lpoly: resid log wage - resid SST, no overlap, more controls
lpoly rwageNO2 rSSTNO2 if tzcode==6 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') noscatter ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") ///
	addplot(lpoly rwageNO2 rSSTNO2 if tzcode==5 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(dash) lcolor(black) || ///
	lpoly rwageNO2 rSSTNO2 if tzcode==2 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(dash_dot) lcolor(black) || ///
	lpoly rwageNO2 rSSTNO2 if tzcode==3 & !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE') & rSSTNO > -.48 & rSSTNO < .41, degree(1) bwidth(`bw2') lpattern(shortdash) lcolor(black)) 
graph export "`work'/graphs/lpoly_rwage_rSST_nooverlap_mc.eps", as(eps) replace
* Binscatter: resid log wage - resid SST, no overlap, more controls
binscatter rwageNO2 rSSTNO2 if !inrange(longitude, `donutPM') & !inrange(longitude, `donutMC') & !inrange(longitude, `donutCE')  & rSSTNO > -.48 & rSSTNO < .41, by(tzcode) colors(black) ///
	legend(lab(1 "Pacific") lab(2 "Mountain") lab(3 "Central") lab(4 "Eastern")) title("") xtitle("Avg. sunset time, residualized") ytitle("ln(Avg. wkly wage), residualized")
graph export "`work'/graphs/binscatter_rwage_rSST_nooverlap_mc.eps", as(eps) replace


* Close graph windows
window manage close graph _all

} /* end flow control for `graphs' */


timer off 1
timer list 1
capture log close



